0) Setup
#libzzs---
library(ggplot2)
library(dplyr)
library(ggpubr)
library(car)
library(emmeans)
library(report)
library(DHARMa)
rm(list=ls())
set.seed(666) #Fix kernell
1) Import data & affect factors
setwd("/Users/martinmartin/Downloads/DATA/Github/Field")
sAll<-read.table("Field_G_A_P_learning.csv", sep="",header=T)
sAll$cat<-as.factor(sAll$cat)
sAll$ID<-as.factor(sAll$ID)
sAll$pck<-as.factor(sAll$pck)
sAll$cond<-as.factor(sAll$cond)
sAll$expe<-as.factor(sAll$expe)
sA1<-sAll %>% filter(cond=="GC"|cond=="PC") %>% mutate(cond="CC")
sA2<-sAll %>% filter(cond!="GC"&cond!="PC")
sAll<-bind_rows(sA1,sA2)
sAll$cond<-as.factor(sAll$cond)
levels(sAll$cond)
## [1] "A-" "A+" "AA" "CC" "G-" "G+" "P-" "P+"
c1<-"A 200\u00b5g/L"
c2<-"A 500\u00b5g/L"
c3<-"Control N°1"
c4<-"Control N°2"
c5<-"G 100\u00b5g/L"
c6<-"G 200\u00b5g/L"
c7<-"P 1mg/L"
c8<-"P 10mg/L"
levels(sAll$cond) <- c(c1,c2,c3,c4,c5,c6,c7,c8)
yo<-sAll %>% count(cond,ID)
sAll$cond<-factor(sAll$cond,levels=c("Control N°1",
"A 200\u00b5g/L",
"A 500\u00b5g/L",
"Control N°2",
"G 100\u00b5g/L",
"G 200\u00b5g/L",
"P 1mg/L","P 10mg/L"))
levels(sAll$cond)
## [1] "Control N°1" "A 200µg/L" "A 500µg/L" "Control N°2" "G 100µg/L"
## [6] "G 200µg/L" "P 1mg/L" "P 10mg/L"
2) Classification
srawAe<-sAll %>%
group_by(cond,expe,ID,cat) %>%
mutate(numbering = row_number()) %>% ungroup()
#Define threshold for being "too much low" to be able to dive
l1<- max(srawAe$pos_y)-max(srawAe$pos_y)/8
#Define poscondon to differentiate individuals that did not respond versus respond
lowl1 = -10
upl1 = +10
scaseAe<-srawAe %>% group_by(cond,ID,cat) %>%
mutate(threshold = case_when(numbering==min(numbering) & pos_y>l1 ~ 1,TRUE~0)) %>% ungroup()
#Count thresholds and sum into "rep" variable
sfilterAe<-scaseAe %>%
group_by(cond,expe,ID,cat) %>%
summarise(rep=sum(threshold,na.rm=T),
dY=sum(-dY,na.rm = T),
pos_y=mean(pos_y,na.rm=T)) %>%
filter(rep==0) %>% ungroup()
#Split positions according to limit
sposAe<-sfilterAe %>%
mutate(slope=case_when(dY < lowl1 ~ "down",
dY>upl1 ~ "up",
dY >= lowl1 & dY <= upl1 ~ "flat"))
s5<-sposAe %>% filter(slope=="up") %>% mutate(rep=1)
s5b<-sposAe %>% filter(slope=="flat") %>% mutate(rep=0)
stotAe<-bind_rows(s5,s5b)
rm(s5,s5b,scaseAe,sposAe,srawAe)
srAe<-stotAe %>%
group_by(cond,cat) %>%
summarise(rep=100*mean(rep,na.rm=TRUE))
stvAe<-sAll %>%
group_by(cond,ID,cat) %>%
summarise(dY=sum(-dY,na.rm = T),
pos_y=mean(pos_y,na.rm=T))
complet<-stvAe %>% group_by(cat) %>% count(cat)
filtered<-stotAe %>% group_by(cat) %>% count(cat)
filteredcond<-stotAe %>% group_by(cond,cat) %>% count(cond,cat)
rm(stvAe,complet,filtered)
3) Count number of individuals –> Supp Table T2
jr<-sAll %>% count(cond,ID)
jcat<-sAll %>% count(cond,ID,cat)
yoae<-jcat %>% group_by(cond,ID) %>% count(ID)
yoae1<-yoae %>% filter(cond==c1)
yoae2<-yoae %>% filter(cond==c2)
yoae3<-yoae %>% filter(cond==c3)
yoae4<-yoae %>% filter(cond==c4)
yoae5<-yoae %>% filter(cond==c5)
yoae6<-yoae %>% filter(cond==c6)
yoae7<-yoae %>% filter(cond==c7)
yoae8<-yoae %>% filter(cond==c8)
aebis<-sfilterAe %>% count(cond,ID)
aebis1<-sfilterAe %>% count(cond,ID,cat)
yoaebis<-aebis1 %>% group_by(cond,ID) %>% count(ID)
aeter<-stotAe %>% count(cond,ID)
aeter1<-stotAe %>% count(cond,ID,cat)
yoaeter<-aeter1 %>% group_by(cond,ID) %>% count(ID)
aeae<-left_join(yoae,yoaebis, by=c("cond"="cond","ID"="ID"))
aeae2<-left_join(yoaebis,yoaeter, by=c("cond"="cond","ID"="ID"))
print(c1)
## [1] "A 200µg/L"
print(length(yoae1$ID))
## [1] 30
aec1<-aeae %>% filter(cond==c1)
print("trial per ID")
## [1] "trial per ID"
sum(aec1$n.x)
## [1] 300
sum(aec1$n.y)
## [1] 234
aec1<-aeae2 %>% filter(cond==c1)
sum(aec1$n.y)
## [1] 234
print(c2)
## [1] "A 500µg/L"
print(length(yoae2$ID))
## [1] 29
aec2<-aeae %>% filter(cond==c2)
print("trial per ID")
## [1] "trial per ID"
sum(aec2$n.x)
## [1] 290
sum(aec2$n.y)
## [1] 225
aec2<-aeae2 %>% filter(cond==c2)
sum(aec2$n.y)
## [1] 224
print(c3)
## [1] "Control N°1"
print(length(yoae3$ID))
## [1] 34
aec3<-aeae %>% filter(cond==c3)
print("trial per ID")
## [1] "trial per ID"
sum(aec3$n.x)
## [1] 340
sum(aec3$n.y,na.rm=T)
## [1] 258
aec3<-aeae2 %>% filter(cond==c3)
sum(aec3$n.y)
## [1] 257
print(c4)
## [1] "Control N°2"
print(length(yoae4$ID))
## [1] 44
aec4<-aeae %>% filter(cond==c4)
print("trial per ID")
## [1] "trial per ID"
sum(aec4$n.x)
## [1] 440
sum(aec4$n.y,na.rm=T)
## [1] 367
aec4<-aeae2 %>% filter(cond==c4)
sum(aec4$n.y)
## [1] 362
print(c5)
## [1] "G 100µg/L"
print(length(yoae5$ID))
## [1] 30
aec5<-aeae %>% filter(cond==c5)
print("trial per ID")
## [1] "trial per ID"
sum(aec5$n.x)
## [1] 300
sum(aec5$n.y,na.rm=T)
## [1] 255
aec5<-aeae2 %>% filter(cond==c5)
sum(aec5$n.y)
## [1] 254
print(c6)
## [1] "G 200µg/L"
print(length(yoae6$ID))
## [1] 30
aec6<-aeae %>% filter(cond==c6)
print("trial per ID")
## [1] "trial per ID"
sum(aec6$n.x)
## [1] 300
sum(aec6$n.y,na.rm=T)
## [1] 216
aec6<-aeae2 %>% filter(cond==c6)
sum(aec6$n.y)
## [1] 214
print(c7)
## [1] "P 1mg/L"
print(length(yoae7$ID))
## [1] 31
aec7<-aeae %>% filter(cond==c7)
print("trial per ID")
## [1] "trial per ID"
sum(aec7$n.x)
## [1] 310
sum(aec7$n.y,na.rm=T)
## [1] 241
aec7<-aeae2 %>% filter(cond==c7)
sum(aec7$n.y)
## [1] 240
print(c8)
## [1] "P 10mg/L"
print(length(yoae8$ID))
## [1] 30
aec8<-aeae %>% filter(cond==c8)
print("trial per ID")
## [1] "trial per ID"
sum(aec8$n.x)
## [1] 300
sum(aec8$n.y,na.rm=T)
## [1] 241
aec8<-aeae2 %>% filter(cond==c8)
sum(aec8$n.y)
## [1] 240
4) Evaluate trial bias –> Supp Table T3
print(c1)
## [1] "A 200µg/L"
R1p<-jcat %>% filter(cond==c1) %>% count(cat)
R1r<-stotAe %>% filter(cond==c1) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
##
## Chi-squared test for given probabilities
##
## data: final1$nz
## X-squared = 13.697, df = 9, p-value = 0.1335
print(c2)
## [1] "A 500µg/L"
R1p<-jcat %>% filter(cond==c2) %>% count(cat)
R1r<-stotAe %>% filter(cond==c2) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
##
## Chi-squared test for given probabilities
##
## data: final1$nz
## X-squared = 11.879, df = 9, p-value = 0.2202
print(c3)
## [1] "Control N°1"
R1p<-jcat %>% filter(cond==c3) %>% count(cat)
R1r<-stotAe %>% filter(cond==c3) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
##
## Chi-squared test for given probabilities
##
## data: final1$nz
## X-squared = 5.3133, df = 9, p-value = 0.8062
print(c4)
## [1] "Control N°2"
R1p<-jcat %>% filter(cond==c4) %>% count(cat)
R1r<-stotAe %>% filter(cond==c4) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
##
## Chi-squared test for given probabilities
##
## data: final1$nz
## X-squared = 11.231, df = 9, p-value = 0.2602
print(c5)
## [1] "G 100µg/L"
R1p<-jcat %>% filter(cond==c5) %>% count(cat)
R1r<-stotAe %>% filter(cond==c5) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
##
## Chi-squared test for given probabilities
##
## data: final1$nz
## X-squared = 13.565, df = 9, p-value = 0.1387
print(c6)
## [1] "G 200µg/L"
R1p<-jcat %>% filter(cond==c6) %>% count(cat)
R1r<-stotAe %>% filter(cond==c6) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
##
## Chi-squared test for given probabilities
##
## data: final1$nz
## X-squared = 7.0233, df = 9, p-value = 0.6347
print(c7)
## [1] "P 1mg/L"
R1p<-jcat %>% filter(cond==c7) %>% count(cat)
R1r<-stotAe %>% filter(cond==c7) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
##
## Chi-squared test for given probabilities
##
## data: final1$nz
## X-squared = 6.2857, df = 9, p-value = 0.711
print(c8)
## [1] "P 10mg/L"
R1p<-jcat %>% filter(cond==c8) %>% count(cat)
R1r<-stotAe %>% filter(cond==c8) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
##
## Chi-squared test for given probabilities
##
## data: final1$nz
## X-squared = 8, df = 9, p-value = 0.5341
5) Plot all trials –> Figure 2) A), B), C)
sf3 <- stotAe %>%
group_by(cond,cat) %>%
summarise(dY=mean(dY,na.rm=T),
rep=mean(rep,na.rm=T)) %>% ungroup()
sf3$cond<-as.factor(sf3$cond)
st1<-sf3 %>% filter(cond=="Control N°1")
st2<-sf3 %>% filter(cond=="A 200\u00b5g/L")
st3<-sf3 %>% filter(cond=="A 500\u00b5g/L")
st4<-sf3 %>% filter(cond=="Control N°2")
st5<-sf3 %>% filter(cond=="P 1mg/L")
st6<-sf3 %>% filter(cond=="P 10mg/L")
st7<-sf3 %>% filter(cond=="G 100\u00b5g/L")
st8<-sf3 %>% filter(cond=="G 200\u00b5g/L")
print("ATRAZINE")
## [1] "ATRAZINE"
ggplot(NULL) +
geom_smooth(data=st1, aes(y=dY, x=cat, group=1, color="Control N°1"), method="gam",
formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
geom_smooth(data=st2, aes(y=dY, x=cat, group=1, color="A 200\u00b5g/L"), method="gam",
formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
geom_smooth(data=st3, aes(y=dY, x=cat, group=1, color="A 500\u00b5g/L"), method="gam",
formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
geom_point(data=st1, aes(y=dY, x=cat, color="Control N°1"), size=3, shape=15) +
geom_point(data=st2, aes(y=dY, x=cat, color="A 200\u00b5g/L"), size=3, shape=15) +
geom_point(data=st3, aes(y=dY, x=cat, color="A 500\u00b5g/L"), size=3, shape=15) +
theme_classic() +
labs(y="Vertical Distance (mm)", x="Trial", color="Condition") +
theme(plot.title = element_text(hjust = 0.5, size=22, face="bold"),
axis.text=element_text(size=20),
axis.title=element_text(size=20)) +
scale_x_discrete(name="Trial",
limits=c("1","2","3","4","5","6","7","8","9","10","Test")) +
scale_color_manual(values = c("Control N°1"="#54C6CC", "A 200\u00b5g/L"="#FFD966",
"A 500\u00b5g/L"="#F0B33E"))
print("PARACETAMOL")
## [1] "PARACETAMOL"
ggplot(NULL) +
geom_smooth(data=st4, aes(y=dY, x=cat, group=1, color="Control N°2"), method="gam",
formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
geom_smooth(data=st5, aes(y=dY, x=cat, group=1, color="P 1mg/L"), method="gam",
formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
geom_smooth(data=st6, aes(y=dY, x=cat, group=1, color="P 10mg/L"), method="gam",
formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
geom_point(data=st4, aes(y=dY, x=cat, color="Control N°2"), size=3, shape=15) +
geom_point(data=st5, aes(y=dY, x=cat, color="P 1mg/L"), size=3, shape=15) +
geom_point(data=st6, aes(y=dY, x=cat, color="P 10mg/L"), size=3, shape=15) +
theme_classic() +
labs(y="Vertical Distance (mm)", x="Trial", color="Condition") +
theme(plot.title = element_text(hjust = 0.5, size=22, face="bold"),
axis.text=element_text(size=20),
axis.title=element_text(size=20)) +
scale_x_discrete(name="Trial",
limits=c("1","2","3","4","5","6","7","8","9","10","Test")) +
scale_color_manual(values = c("Control N°2"="#54C6CC", "P 1mg/L"="#00A1EA",
"P 10mg/L"="#0070C0"))
print("GLYPHOSATE")
## [1] "GLYPHOSATE"
ggplot(NULL) +
geom_smooth(data=st4, aes(y=dY, x=cat, group=1, color="Control N°2"), method="gam",
formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
geom_smooth(data=st7, aes(y=dY, x=cat, group=1, color="G 100\u00b5g/L"), method="gam",
formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
geom_smooth(data=st8, aes(y=dY, x=cat, group=1, color="G 200\u00b5g/L"), method="gam",
formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
geom_point(data=st4, aes(y=dY, x=cat, color="Control N°2"), size=3, shape=15) +
geom_point(data=st7, aes(y=dY, x=cat, color="G 100\u00b5g/L"), size=3, shape=15) +
geom_point(data=st8, aes(y=dY, x=cat, color="G 200\u00b5g/L"), size=3, shape=15) +
theme_classic() +
labs(y="Vertical Distance (mm)", x="Trial", color="Condition") +
theme(plot.title = element_text(hjust = 0.5, size=22, face="bold"),
axis.text=element_text(size=20),
axis.title=element_text(size=20)) +
scale_x_discrete(name="Trial",
limits=c("1","2","3","4","5","6","7","8","9","10","Test")) +
scale_color_manual(values = c("Control N°2"="#54C6CC", "G 100\u00b5g/L"="#EB8176",
"G 200\u00b5g/L"="#941100"))
6) Plot trials 1 and Test –> Figure 2 D)
sx<-stotAe %>%
filter(cat==1|cat==10) %>%
group_by(cond,ID,cat) %>%
summarise(dY=mean(dY,na.rm = T),
rep=mean(rep,na.rm=T)) %>% ungroup()
sx %>% group_by(cond,cat) %>% count()
## # A tibble: 16 × 3
## # Groups: cond, cat [16]
## cond cat n
## <fct> <fct> <int>
## 1 Control N°1 1 29
## 2 Control N°1 10 26
## 3 A 200µg/L 1 29
## 4 A 200µg/L 10 26
## 5 A 500µg/L 1 23
## 6 A 500µg/L 10 20
## 7 Control N°2 1 40
## 8 Control N°2 10 31
## 9 G 100µg/L 1 28
## 10 G 100µg/L 10 21
## 11 G 200µg/L 1 24
## 12 G 200µg/L 10 22
## 13 P 1mg/L 1 27
## 14 P 1mg/L 10 25
## 15 P 10mg/L 1 27
## 16 P 10mg/L 10 21
dae<-glm(dY~cat*cond, data=sx)
dx<-sx %>% count(cond,cat)
yv <- predict(dae,type="link",se.fit=TRUE)
ydAE<-data.frame(yv$fit,yv$se.fit,sx$cat,sx$cond)
yd2<-ydAE %>% rename(cond=sx.cond) %>% rename(cat=sx.cat)
inv<-family(dae)$linkinv
ggplot(sx, aes(x=cat,y=dY))+
geom_pointrange(data=yd2,aes(x=cat,y=inv(yv.fit),
ymin=inv(yv.fit-1.96*yv.se.fit),
ymax=inv((yv.fit+1.96*yv.se.fit)),
group=cond,color=cond),size=2,linewidth=3)+
geom_point(aes(color=cond))+
theme_classic() +
labs(y="Vertical Distance (mm)",x="Trial")+
theme(plot.title = element_text(hjust = 0.5))+
theme(plot.title = element_text(size=22,face = "bold"))+
theme(axis.text=element_text(size=20,color="black"),
axis.title=element_text(size=20,color="black"))+
scale_color_manual(values=c("#54C6CC","#FFD966","#F0B33E","#54C6CC","#EB8176","#941100",
"#00A1EA","#0070C0"))+
theme(legend.title = element_text(size=20,color="black"),
legend.text = element_text(size=20))+
guides(colour=guide_legend(title="Species"))+
theme(legend.position = "none",
strip.text = element_text(size = 18))+
stat_compare_means(comparisons = list(c("1","10")),
aes(label = ..p.signif..),
bracket.size = 1,size=5)+
facet_wrap(~cond,nrow=1)
7) Stats
SR1<-sx %>% filter(cond==c1)
mR1<-lmerTest::lmer(dY~cat+(1|ID),data=SR1)
simR1 <- simulateResiduals(fittedModel = mR1, plot = T)
Anova(mR1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dY
## Chisq Df Pr(>Chisq)
## cat 6.9706 1 0.008286 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR1, pairwise ~ cat,adjust="tukey")
## $emmeans
## cat emmean SE df lower.CL upper.CL
## 1 29.1 2.23 53 24.6 33.5
## 10 20.5 2.37 53 15.7 25.2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## cat1 - cat10 8.57 3.26 27.7 2.631 0.0137
##
## Degrees-of-freedom method: kenward-roger
SR2<-sx %>% filter(cond==c2)
mR2<-lmerTest::lmer(dY~cat+(1|ID),data=SR2)
simR2 <- simulateResiduals(fittedModel = mR2, plot = T)
Anova(mR2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dY
## Chisq Df Pr(>Chisq)
## cat 4.0355 1 0.04455 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR2, pairwise ~ cat,adjust="tukey")
## $emmeans
## cat emmean SE df lower.CL upper.CL
## 1 31.5 2.67 40.9 26.1 36.9
## 10 24.0 2.88 41.0 18.2 29.8
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## cat1 - cat10 7.46 3.78 23.2 1.974 0.0604
##
## Degrees-of-freedom method: kenward-roger
SR3<-sx %>% filter(cond==c3)
mR3<-lmerTest::lmer(dY~cat+(1|ID),data=SR3)
simR3 <- simulateResiduals(fittedModel = mR3, plot = T)
Anova(mR3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dY
## Chisq Df Pr(>Chisq)
## cat 11.613 1 0.0006549 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR3, pairwise ~ cat,adjust="tukey")
## $emmeans
## cat emmean SE df lower.CL upper.CL
## 1 29.7 2.03 49.9 25.6 33.8
## 10 21.4 2.14 51.3 17.1 25.7
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## cat1 - cat10 8.28 2.46 26.7 3.373 0.0023
##
## Degrees-of-freedom method: kenward-roger
SR4<-sx %>% filter(cond==c4)
mR4<-lmerTest::lmer(dY~cat+(1|ID),data=SR4)
simR4 <- simulateResiduals(fittedModel = mR4, plot = T)
Anova(mR4)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dY
## Chisq Df Pr(>Chisq)
## cat 7.9536 1 0.004799 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR4, pairwise ~ cat,adjust="tukey")
## $emmeans
## cat emmean SE df lower.CL upper.CL
## 1 26.5 1.85 69 22.9 30.2
## 10 18.7 2.12 69 14.4 22.9
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## cat1 - cat10 7.88 2.81 38.5 2.801 0.0079
##
## Degrees-of-freedom method: kenward-roger
SR5<-sx %>% filter(cond==c5)
mR5<-lmerTest::lmer(dY~cat+(1|ID),data=SR5)
simR5 <- simulateResiduals(fittedModel = mR5, plot = T)
Anova(mR5)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dY
## Chisq Df Pr(>Chisq)
## cat 1.1381 1 0.2861
emmeans(mR5, pairwise ~ cat,adjust="tukey")
## $emmeans
## cat emmean SE df lower.CL upper.CL
## 1 22.7 2.54 47 17.6 27.8
## 10 18.5 2.97 47 12.6 24.5
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## cat1 - cat10 4.14 3.91 25.2 1.059 0.2995
##
## Degrees-of-freedom method: kenward-roger
SR6<-sx %>% filter(cond==c5)
mR6<-lmerTest::lmer(dY~cat+(1|ID),data=SR6)
simR6 <- simulateResiduals(fittedModel = mR6, plot = T)
Anova(mR6)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dY
## Chisq Df Pr(>Chisq)
## cat 1.1381 1 0.2861
emmeans(mR6, pairwise ~ cat,adjust="tukey")
## $emmeans
## cat emmean SE df lower.CL upper.CL
## 1 22.7 2.54 47 17.6 27.8
## 10 18.5 2.97 47 12.6 24.5
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## cat1 - cat10 4.14 3.91 25.2 1.059 0.2995
##
## Degrees-of-freedom method: kenward-roger
SR7<-sx %>% filter(cond==c7)
mR7<-lmerTest::lmer(dY~cat+(1|ID),data=SR7)
simR7 <- simulateResiduals(fittedModel = mR7, plot = T)
Anova(mR7)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dY
## Chisq Df Pr(>Chisq)
## cat 6.9368 1 0.008444 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR7, pairwise ~ cat,adjust="tukey")
## $emmeans
## cat emmean SE df lower.CL upper.CL
## 1 25.0 2.40 45.7 20.1 29.8
## 10 17.6 2.49 46.9 12.6 22.6
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## cat1 - cat10 7.37 2.81 25.1 2.618 0.0148
##
## Degrees-of-freedom method: kenward-roger
SR8<-sx %>% filter(cond==c8)
mR8<-lmerTest::lmer(dY~cat+(1|ID),data=SR8)
simR8 <- simulateResiduals(fittedModel = mR8, plot = T)
Anova(mR8)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dY
## Chisq Df Pr(>Chisq)
## cat 2.0166 1 0.1556
emmeans(mR8, pairwise ~ cat,adjust="tukey")
## $emmeans
## cat emmean SE df lower.CL upper.CL
## 1 26.4 2.51 37.6 21.3 31.4
## 10 22.6 2.76 42.7 17.1 28.2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## cat1 - cat10 3.73 2.65 21.8 1.408 0.1732
##
## Degrees-of-freedom method: kenward-roger